Assessing spatial distribution, genetic variants, and virulence of pathogen Mycoplasma agassizii in threatened Mojave desert tortoises

Abstract Mojave desert tortoises (Gopherus agassizii), a threatened species under the US Endangered Species Act, are long‐lived reptiles that experience a chronic respiratory disease. The virulence of primary etiologic agent, Mycoplasma agassizii, remains poorly understood, but it exhibits temporal and geographic variability in causing disease outbreaks in host tortoises. Multiple attempts to culture and characterize the diversity of M. agassizii have had minimal success, even though this opportunistic pathogen chronically persists in nearly every population of Mojave desert tortoises. The current geographic range and the molecular mechanisms of virulence of the type‐strain, PS6T, are unknown, and the bacterium is thought to have low‐to‐moderate virulence. We designed a quantitative polymerase chain reaction (qPCR) targeting three putative virulence genes annotated on the PS6T genome as exo‐α‐sialidases, enzymes which facilitate growth in many bacterial pathogens. We tested 140 M. agassizii‐positive DNA samples collected from 2010 to 2012 across the range of Mojave desert tortoises. We found evidence of multiple‐strain infections within hosts. We also found the prevalence of these sialidase‐encoding genes to be highest in tortoise populations surrounding southern Nevada, the area from which PS6T was originally isolated. We found a general pattern of loss or reduced presence of sialidase among strains, even within a single host. However, in samples that were positive for any of the putative sialidase genes, one particular gene (528), was positively associated with bacterial loads of M. agassizii and may act as a growth factor for the bacterium. Our results suggest three evolutionary patterns: (1) high levels of variation, possibly due to neutral changes and chronic persistence, (2) a trade‐off between moderate virulence and transmission, and (3) selection against virulence in environmental conditions known to be physiologically stressful to the host. Our approach of quantifying genetic variation via qPCR represents a useful model of studying host–pathogen dynamics.


| INTRODUC TI ON
Amid increasing urgency to develop standards in monitoring pathogen variation and movement across landscapes, wildlife biologists play an important role in developing pathogen-specific methodology to understand ongoing threats of emerging diseases (Watsa and WDSFG, 2020). Many standard methods for surveillance of pathogen variation, including next-generation sequencing (NGS), whole-genome shotgun (WGS) sequencing, and sequencing of cultured isolates, are often biased in detecting the dominant or fastestgrowing strains while missing rare species or strains (Dyrhovden et al., 2019;Varghese et al., 2017). Alternatively, probe-based quantitative PCR (qPCR) can be very specific, which makes it ideal for targeted pathogen surveys even if its specificity can result in false negatives due to slight changes in the target sequence (Braun et al., 2014;Geraci-Yee et al., 2022;Rockett et al., 2011).
For Mycoplasma agassizii, the pathogen that causes upper respiratory tract disease (URTD) in Mojave desert tortoises (Gopherus agassizii; Brown et al., 1994), culture methods and antibody testing have had limited effectiveness in diagnostic testing and identifying strain variation (Braun et al., 2014;Luzuriaga-Neira et al., 2021;Sandmeier, Weitzman, et al., 2017). In thirty years since its first isolation from its desert tortoise hosts, only two strains of M. agassizii have been isolated from this tortoise species (Brown et al., 1994(Brown et al., , 2001Jacobson et al., 1991). Hydrolysis probe-based quantitative PCR has been instrumental in diagnosing M. agassizii in tortoises, as loads of M. agassizii DNA commonly are too low to sequence accurately (Braun et al., 2014). These assays have also led to insights on ecological patterns of disease (Aiello et al., 2014;Braun et al., 2014;Sandmeier et al., 2018). For example, low loads of M. agassizii often do not cause disease, and M. agassizii persists in tortoises for many years (Jacobson et al., 2014;Sandmeier et al., 2018;Weitzman et al., 2017). Recrudescence of M. agassizii loads within hosts and transmission from other tortoises both are known to cause increases in disease prevalence in populations (Aiello et al., 2016;. However, the mechanisms causing this variation in growth is unknown and no clear virulence factors have been identified for this bacterium (Luzuriaga-Neira et al., 2021;Sandmeier et al., 2018).
Virulence of M. agassizii is thought to be relatively low in desert tortoises in years in which environmental conditions do not deviate from mean conditions, outbreaks of disease are relatively rare, and levels of morbidity are low (Jacobson et al., 2014;Sandmeier et al., 2009;USFWS, 2015). Virulence is associated with adaptive antibody responses in the host, which appear to maintain local inflammation of the mucosal membranes . Inflammation and nasal exudate seem to increase the chances of transmission among hosts, linking virulence and transmission (Alizon et al., 2008;Ebert & Bull, 2008). Importantly, innate immune mechanisms appear to be more effective at reducing loads of M. agassizii than adaptive immune mechanisms, which is common in reptiles Zimmerman et al., 2010). Success of new variants of M. agassizii likely is a balance between transmission and virulence, mediated by these host immune responses (Sandmeier & Tracy, 2014).
Recent genetic analysis of M. agassizii from wild tortoises found low amounts of genetic variation and confirmed that the typestrain from the Mojave Desert, PS6 T (ATCC® 700616 T ), appears to have persisted in desert tortoise populations after its initial isolation in the early 1990s (Brown et al., 1994(Brown et al., , 2001Luzuriaga-Neira et al., 2021;Weitzman et al., 2017). To understand virulence and genetic variation of PS6 T , we identified virulence gene candidates to target with qPCR using data obtained by the recent sequencing and annotation of this strain (Alvarez-Ponce et al., 2018). We targeted three annotated exoα-sialidase genes found in PS6 T , because sialidases are enzymes identified as important virulence factors in similar bacterial pathogens (Brown et al., 2004;Severi et al., 2007).
Bacterial sialidases can function as virulence factors by freeing host monosaccharides called sialic acids, abundant throughout the mucosa of deuterostomes (Corfield, 1992;Lewis & Lewis, 2012). In doing so, they make additional nutritional sources available for bacterial growth, expose new sites for adhesion on host substrate, and increase biofilm formation (Hardy et al., 2017;Severi et al., 2007;Vimr et al., 2004). Therefore, we predicted that sialidases of M. agassizii increase bacterial growth rates and lead to increased disease (Michaels et al., 2019;Sandmeier et al., 2018;Weitzman et al., 2017).
The Mojave desert tortoise-M. agassizii system provides an opportunity to analyze low-to-moderate virulence under stabilized conditions, where a chronic, opportunistic pathogen evolved with slow immune responses and slow transmission rates of its host (Aiello et al., 2016;Luzuriaga-Neira et al., 2021;Sandmeier et al., 2018). Table 1 outlines multiple, non-mutually exclusive hypotheses of how virulence genes may evolve and affect tortoises and ecological patterns across the landscape (Dykhuizen & Kalia, 2008;Ebert & Bull, 2008). We made several predictions about how these patterns may manifest in desert tortoise populations (Table 1). We expected a K E Y W O R D S exoα-sialidase, host-pathogen dynamics, mycoplasma, testudines, virulence gene, wildlife disease

T A X O N O M Y C L A S S I F I C A T I O N
Applied ecology, Conservation ecology, Conservation genetics, Disease ecology, Ecological genetics, Evolutionary ecology, Functional ecology, Genetics, Microbial ecology, Spatial ecology TA B L E 1 Theoretical patterns of evolution of virulence, specific to the Mojave desert tortoise-Mycoplasma agassizii system.  Dykhuizen and Kalia (2008) and  5. Putative virulence genes have evolved to avirulence Neither presence nor loads of the genes are associated with URTD

No regional pattern
Unlikely: sialidase genes are often associated with virulence in mycoplasmas and other opportunistic pathogens Brown et al. (2001) and Severi et al. (2007) Note: We detail expectations for each pattern in individuals across the geographic landscape, and whether we deem these patterns to likely occur in our host-pathogen system. trade-off between virulence and transmission rates, virulence to be affected by host physiology and immunocompetence, and possible multiple-strain infections ( Table 1).
We used samples of M. agassizii collected via lavages the nasal cavities of wild tortoises to quantify genetic gains and losses of three putative exo-sialidase genes, each of which is represented by one unique copy on the PS6 T genome (Alvarez-Ponce et al., 2018). We used qPCR of the single-copy 16S rRNA gene on the genome of all strains of M. agassizii to quantify total amounts of M. agassizii within samples (Braun et al., 2014;Weitzman et al., 2017). Samples were collected over a three-year period (2010-2012) for health assessment of tortoise populations across the majority of the geographic extent of the Mojave Desert ( Figure 1; described in Sandmeier et al., 2018;Weitzman et al., 2017). Throughout this article, we used "strain" interchangeably with genetic variant, to indicate quantitative differences of these virulence genes compared to PS6 T (Balmer & Tanner, 2011).
First, we confirmed that PS6 T sialidase genes were annotated correctly and likely to be functional by comparing them to genetic and amino acid sequences of known sialidases (Roggentin et al., 1993;Taylor, 1996;Vimr, 1994). We assessed how presence and loads of the three genes were associated and then tested four hypotheses.
(1) We hypothesized that we would be able to detect strain variation of M. agassizii across the range of Mojave desert tortoises and within samples. If we found a partial loss of these genes within samples, this would indicate multiple-strain infections in hosts. (2) We tested whether the presence of sialidase genes was associated with three measures of URTD: mycoplasmal load, signs of disease (nasal exudate), and disease severity Sandmeier et al., 2018).
We predicted that detection of the genes, as well as the number of the three genes detected in each sample (0-3), would predict greater M. agassizii loads, presence of disease, and disease severity (Roggentin et al., 1993;Taylor, 1996;Vimr, 1994). Because sialidase genes clustered on the landscape and were absent from many samples (Figure 1), we also analyzed the subsets of samples in which at least one putative sialidase gene was present.
(3) We hypothesized that we would find a relationship between genetic loads of the sialidase genes when they were present and the three measures of URTD (M. agassizii load, disease presence, and disease severity). (4) Lastly, since Mojave desert tortoises are separated into three distinct genotypic populations driven by geographic distance, we hypothesized that prevalence of the sialidase genes would vary among these distinct host populations (Hagerty & Tracy, 2010;Shaffer et al., 2017).

| Study population and data collection
Across the entire Mojave Desert, 20 populations of wild tortoises (n = 419) were sampled during their active seasons (April-June) in 2010-2012 (described in Sandmeier et al., 2018;Weitzman et al., 2017). Individual tortoises were evaluated for presence (0 or 1) and severity (scored 0-6) of signs of disease Weitzman et al., 2017). Nasal lavage samples were collected from each tortoise by rinsing 3 mL of sterile saline solution through the nares. The lavage samples were preserved in RNAlater stabilizing solution (Ambion Inc.) and placed on ice until they could be frozen within 12 h of collection . DNA was then extracted from 500 μL of the samples with the Qiagen DNeasy Blood and Tissue Kit (Qiagen Inc.) and stored at −20°C . One-hundred forty samples from 17 local populations tested positive for M. agassizii by qPCR Weitzman et al., 2017) and had enough volume of DNA sample remaining to be tested by four separate additional assays: three assays for each of the putative sialidase genes, and one assay for the quantification of M. agassizii (Braun et al., 2014).
We conducted BLASTp searches in GenBank to examine the accuracy of the annotation of the three genes by comparing sequence similarity to sialidase genes known to be functional in other bacteria. Using T-Coffee (Di Tommaso et al., 2011), we aligned amino acid sequences of each M. agassizii PS6 T sialidase gene with other F I G U R E 1 Distribution of Mycoplasma agassizii-positive tortoise samples (n = 140) used in this study. Quantitative PCR identified 31 samples with at least one sialidase gene (black points), mainly clustered in southern Nevada, where PS6 T was originally isolated (Brown et al., 1994(Brown et al., , 2001. known bacterial sialidases to verify presence of conserved amino acid motifs, found in all functional bacterial sialidases: the YRIP motif (Arg-Ile/Leu-Pro) and the Asp-box motif (Ser/Thr-X-Asp-[X]-Gly-X-Thr-Trp/Phe) (Appendix S1; Roggentin et al., 1993;Taylor, 1996;Vimr, 1994).

| Quantitative PCR assay design and optimization
Primers and probe candidates were generated with Primer Express Software v3.0.1 under default parameters for ideal qPCR chemistry (Thermo Fisher Scientific™). Target amplicons were verified to target CDS regions using NCBI GenBank, and BLASTn queries verified specificity for each Primer Express output. Specificity analyses using BLASTn eliminated any amplicon candidates with high similarity to non-target sequences of other microbes. Each final amplicon with primer and probe sequences was selected to be specific to each PS6 T sialidase gene (  700618)) and the other target sialidase sequences. We also tested assays against one plate of M. agassiziinegative samples that were collected from the same initial disease survey Weitzman et al., 2017). All Ct values were converted to the number of copies of the target DNA. We used synthetic gBlock gene fragments as ampliconspecific positive controls ( Table 2; Integrated DNA Technologies, San Diego, California, USA). The gBlock positive controls were run on each plate in triplicate in 10-fold serial dilution curves, over a concentration range of 5 × 10 6 copies to 5 copies along with negative controls.

| Statistical analyses
We first tested whether the presence (0 = absent, 1 = present) of each of the three genes was correlated with each other via pairwise chi-square analyses. Within the samples positive for any sialidase gene, we also used a correlation matrix to test for associations among loads of the three genes. We then tested for a relationship between total number of sialidase genes per sample (0-3) and mean load of    query cover, 48% identities). When the two genes were aligned using BLASTp, the results show even higher similarity (E-value 8 × 10 −159 , 96% query cover, 46% identities; Appendix S1).

| Sialidase genes qPCR
We located conserved amino acid motifs which are present on all functional bacterial sialidases: all three genes contain the YRIP motif, which occurs in the active site of bacterial sialidases (Roggentin et al., 1993). Genes 528 and 905 have one repeat of the Asp-box motif, which likely functions in protein folding for the β-propellor domain, an important structural element in bacterial sialidases (Appendix S1; Quistgaard & Thirup, 2009;Roggentin et al., 1993;Vimr, 1994). We detected genetic variations in the Asp-box motif of gene 906. It is unclear whether these variations represent loss-of-function mutations, but because all functional bacterial sialidases contain both motifs, we were unable to verify likely functionality of gene 906 (Appendix S1; Quistgaard & Thirup, 2009;Roggentin et al., 1993;Vimr, 1994).
In silico and in vitro specificity tests verified that the qPCR assays amplified only the target genes from M. agassizii. The specific probe-primer  Table 2 reports average parameters for each of the three qPCR assays based on standard curves.

| Statistical analyses
Out of these 140 M. agassizii-positive samples, 28 were positive for gene 528, 14 were positive for gene 905, and eight were positive for gene 906 (Figures 1 and 2). Samples tended to have more than one sialidase gene, and we found that the presence of these genes correlated with each other significantly (Figure 2a Finally, loads of the three genes were not associated with disease severity (528: partial p = .627; 905: partial p = .150; 906: partial p = .110).
Total load of sialidase genes significantly varied among the three genetic populations of tortoise hosts (F 2,137 = 4.800, p = .010; Figure 5). There was no difference between the "Northeast Mojave" and the other two genetic populations (Las Vegas: p = .268; California: p = .438), but there was a significant higher level of sialidase genes in the "Las Vegas" versus the "California" populations (p = .007; Figure 5).
Across the geographic range, sialidase genes were especially prevalent around the Las Vegas area and generally absent from southwestern California and the northeastern Mojave (Figures 1 and 3).

| Summary of genetic variation within and across tortoise hosts
Specific to our hypotheses, we found that a variety of M. agassizii genetic strains occur across the Mojave as indicated by the F I G U R E 3 Prevalence of disease and prevalence of sialidase genes among population groups of Mojave desert tortoises (n ≥ 5 for labeled groups) within three distinct genetic populations (Hagerty & Tracy, 2010). Pie chart locations represent the approximate centroids of the sampled individuals in a valley. Visible clusters of sialidase genes do not always appear coincident with URTD. Solid color in the pie charts indicates the proportion positive for either measure. absence of PS6 T sialidase genes in a majority of samples (Figure 2b).
Similarly, we detected multiple-strain infections within individual hosts (Figure 2b). We found mixed support for our hypotheses of a relationship between sialidase genes and virulence based on three measures of disease: mycoplasmal load, presence of disease, and severity of signs of disease. Low R 2 values suggest that the biological relevance of observed relationships is weak. However, in sialidase-positive samples, load of gene 528 was strongly associated with higher loads of M. agassizii and loads of all three genes were strongly associated with increased presence of URTD, suggesting important within-host dynamics involving growth rates and other mechanisms of virulence such as evasion of immune responses (Table 1). Spatially, we found an absence of these sialidase genes in the southwestern and northeastern extremes of the Mojave Desert and a distribution centered around southern Nevada (Figures 1 and 3). To summarize, the type-strain PS6 T is not highly virulent , but higher loads of all three putative sialidase genes are associated with an increased likelihood of mild signs of disease when they are present Weitzman et al., 2017).
Most samples with high loads of sialidase genes contained all three genes, suggesting a synergism in functionality (Figure 2a).
Both genes 528 and 905 are likely functional sialidase genes based on in silico comparisons to known sialidase genes in other organisms. Of the three targeted PS6 T sialidase genes, 528 was present in the most samples and seems to play a role in bacterial nutrition, as its load was linked to higher loads of M. agassizii (Figure 4a), a predictor of recurrent URTD in tortoises and higher prevalence of disease in populations (Aiello et al., 2016;. Gene 528 is also a homolog of the M. alligatoris nanI protein, which increases mycoplasmal growth in highly virulent mycoplasmal infections in alligators (Michaels et al., 2019). Genes 905 and 906 were not associated with increased loads of M. agassizii, but both were associated with disease, suggesting other functions such as providing potential sites for adhesion on degraded mucosal substrate or avoiding immune detection (Corfield, 1992;Taylor, 1996;Vimr et al., 2004). When sialidase activity frees host molecules, this can allow bacteria to mimic host tissue and evade immune response by incorporating host sugars into their own membrane (Vimr, 1994).
Gene 906 may lack sialidase functionality based on in silico analyses and the current understanding of bacterial sialidase gene structure, and we do not view it as a functional gene (Appendix S1; Roggentin et al., 1993;Taylor, 1996;Vimr, 1994). In addition, it occurred only in eight samples that also contained genes 528 or 905, complicating the interpretation of its effect on disease.

| Multiple-strain infections
The multiple-strain infections that we documented here could have occurred by direct transmission of multiple strains or from multiple transmission events, as M. agassizii is known to persist for years in healthy Mojave desert tortoises Sandmeier, Weitzman, et al., 2017). Alternatively, within-host evolution could also lead to a diversity of strains (Table 1; Ebert & Bull, 2008;Sandmeier & Tracy, 2014). High levels of genetic diversity due to genetic drift, within and among individual hosts, can be a hallmark of chronic opportunistic pathogens which generally co-exist with their hosts (Ebert & Bull, 2008 Figure 2b; Sandmeier et al., 2018).
There may be important ecological implications of multiple-strain infections, largely based on effects of microbial interactions within the host, although effects on virulence are often unpredictable in such systems (Balmer & Tanner, 2011;Sofonea et al., 2017). The fact that many tortoises had lower levels of sialidase genes than M. agassizii suggests that even a small number of strains in a tortoise's respiratory tract can free nutrients, potentially benefiting other strains of M. agassizii or other bacterial species in the same location (Lewis & Lewis, 2012;Vimr, 1994).

| Virulence-transmission balance
We interpret our results as a possible balance between M. agassizii persistence in a host at low enough levels to evade innate immune responses and higher growth rates that cause inflammatory responses and higher transmission events ( Table 1; Aiello et al., 2016).
An interesting aspect of this host-pathogen system is that M. agassizii and URTD have been shown to be negatively associated with innate immunity-specifically, high levels of B1 lymphocytes which produce natural antibodies and phagocytize bacterial cells (Sandmeier et al., 2013Slama et al., 2022). Low levels of M. agassizii likely persist in the host without innate response completely clearing the infection Slama et al., 2022). Conversely, high levels of inflammation and nasal exudate seem only temporarily effective in reducing URTD in desert tortoises and seem to maintain chronic disease and risk of transmission (Aiello et al., 2016;. The low levels of sialidase genes in many samples (Figure 2b) may represent a balance of virulence and transmission that primarily results in low-level morbidity, but increased transmission rates. The duplication of sialidase genes in a handful of samples may imply the importance of these genes in transmission, as mycoplasmas have minimal genomes in which gene loss is more common than gene duplication (Razin, 1996;Razin & Hayflick, 2009).
Virulence genes may not exist in areas in which hosts experience conditions known to negatively impact tortoise physiology, such as drier and colder weather regimes, due to past host mortality events (  Figure 3). For example, tortoises in southern California historically experienced the most extreme drought events, the largest decreases in population densities, and the only well-documented disease outbreaks associated with high mortality rates (USFWS, 2015).
A highly virulent strain of M. agassizii was isolated from one population, "Fremont Kramer" (Figure 3), during a disease outbreak (Brown et al., 1994;Jacobson et al., 1991), and only low prevalence of the sialidase genes was detected in the same area during the 2010-2012 survey. Tortoises in the northernmost population in Utah experience a colder climate and hibernate communally in caves, factors thought to promote transmission of M. agassizii and disease (Sandmeier et al., 2013, and we detected no sialidase genes in this area (Figures 1 and 3). In both cases, more virulent strains may overwhelm individual tortoises, tipping virulence levels to higher host mortality rates, thus leading to local extirpation of virulent strains ( Table 1). comm.). Samples positive for M. agassizii could be tested for sialidase genes to confirm localized impacts of loads of these genes on disease (Figures 3 and 4). Assays for sialidase genes 528 and 905 could provide a tool for regional disease and health management in tortoise populations in regions of southern NV regions and eastern CA ( Figure 3), especially in times of increased environmental stochasticity and sudden regional changes in disease prevalence (Averill-Murray et al., 2021;Berry et al., 2020;Peterson, 1996).
This study demonstrated the strengths of using qPCR to study intraspecific microbial diversity in uncultured samples by targeting putative virulence genes and showed that qPCR can provide a sensitive technique to quantify gene-specific variation across a host species' range (Brown et al., 2001;Hardy et al., 2017;Luzuriaga-Neira et al., 2021). However, it is highly likely that our gene-specific qPCRs did not detect other sialidase genes, as even small mutations affect detection and accuracy of highly specific qPCR assays (Braun et al., 2014;Rockett et al., 2011;Vanysacker et al., 2014). Developing similar assays for additional putative virulence genes of PS6 T , and for other strains as they are sequenced, may provide an increasingly more complete picture of varying levels of disease risk across the landscape.
We suspect that difficulties culturing M. agassizii may be influenced by presence of sialidase. The two M. agassizii strains successfully cultured in the past, PS6 T (Brown et al., 2001;Jacobson et al., 1991) and an unnamed California strain no longer available to researchers (Brown et al., 1994), were grown in SP-4 medium, which is rich in sialic acids due to compounds such as fetal bovine serum (Brown et al., 2001;Corfield, 1992;Jacobson et al., 1991;Sherblom et al., 1988). These previous studies tested several types of media, and growth success was only achieved using SP-4 (Brown et al., 2001;Jacobson et al., 1991). Successful culturing of unknown strains which are clearly present in tortoises may require bound sialic acids within media to be broken down, as some strains may have the ability to catabolize sialic acid even without their own sialidase gene (Lewis & Lewis, 2012;Vimr, 1994). Perhaps supplementing sialidase enzyme in order to free up sialic acids bound in the culture medium may allow isolation of a greater diversity of strains.

| CON CLUS ION
A greater diversity of pathogen-surveillance methods will increase knowledge of the mechanisms causing disease, both in threatened species such as Mojave desert tortoises, and other taxa of wildlife (Watsa and WDSFG, 2020). We demonstrated that qPCR methods may be adapted for tracking genetic variants, which may be a more accessible method to survey targeted genetic variation across a large geographic expanse and especially effective for identifying variation in bacteria which occur in multiple-strain infections or are difficult to culture. We were also able to use this technique to suggest hypotheses for the evolution and current prevalence of virulent and avirulent strains. Our study design can easily be adapted for use in other host-pathogen systems to better understand functional dynamics, especially of chronic diseases with complex etiologies.

ACK N OWLED G M ENTS
The authors thank the US Fish and Wildlife Service (Grant F19AC00070) for primary funding of this research. We are also grateful to Colorado State University-Pueblo and the University's CBASE program for providing equipment and software for laboratory work as well as funding undergraduate student salaries (Communities to Build Active STEM Engagement, Department of Education Title III, Award F19AC00070).

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are provided for peer review and are shared publicly in a repository: https://doi.org/10.6084/m9.figsh are.20520 771.v1.